clc
clear

stencil_width = 5;

syms dx
q = sym('q',[stencil_width 1]);

qL = WENO5(q,-1,dx);
qR = WENO5(q, 1,dx);

for iCell = 1:stencil_width
    dqLdq(iCell) = diff(qL,q(iCell));
    dqRdq(iCell) = diff(qR,q(iCell));
end

ccode([dqLdq;dqRdq],'file','WENO5_jab.c');

% for iCell = 1:stencil_width
%     for jCell = 1:stencil_width
%         d2qLdq2(iCell,jCell) = diff(dqLdq(iCell),q(jCell));
%         d2qRdq2(iCell,jCell) = diff(dqRdq(iCell),q(jCell));
%     end
% end
% ccode([d2qLdq2;d2qRdq2],'file','WENO5_jab2.c');

function qrec = WENO5(Q,dir,dx)
syms epsilon
% c = sym('C',[1 3]);
c(1:3,1) = [0.1,0.6,0.3];

stenil_coef (1,:) = [1./3., -7./6., 11./6.,    0.,     0.];
stenil_coef (2,:) = [   0., -1./6.,  5./6., 1./3.,     0.];
stenil_coef (3,:) = [   0.,     0.,  1./3., 5./6., -1./6.];
stenil_coefA(1,:) = [  1., -2.,  1.,  0.,  0.];
stenil_coefA(2,:) = [  0.,  1., -2.,  1.,  0.];
stenil_coefA(3,:) = [  0.,  0.,  1., -2.,  1.];
stenil_coefB(1,:) = [  1., -4.,  3.,  0.,  0.];
stenil_coefB(2,:) = [  0., -1.,  0.,  1.,  0.];
stenil_coefB(3,:) = [  0.,  0.,  3., -4.,  1.];
        
if(dir>0)
    Qr = Q;
else
    Qr = Q(5:-1:1);
end

stencil = stenil_coef  * Qr;
coefA   = stenil_coefA * Qr;
coefB   = stenil_coefB * Qr;

beta = coefA.^2 * 13. / 12. + coefB.^2 * 0.25;

% omega = WENO_JS_weights(c,beta,epsilon);
omega = WENO_Z_weights(c,beta,epsilon);
% omega = WENO_M_weights(c,beta,epsilon);

qrec = sum(stencil.*omega);

end

function w = WENO_Z_weights(c,beta,epsilon)

tau  = abs( beta(3) - beta(1) );
xi = ( tau + epsilon ) ./ ( beta + epsilon );

alpha = c .* ( 1. + xi.^2 );

sum_alpha = sum(alpha);

w = alpha ./ sum_alpha;
end

function w = WENO_M_weights(c,beta,eps)

w_js = WENO_JS_weights(c,beta,eps);

g = w_js .* ( c + c.^2 - 3 .* c .* w_js + w_js.^2 ) ./ ( c.^2 + w_js .* ( 1. - 2. * c ) );

sum_g = sum(g);

w = g ./ sum_g;

end

function w = WENO_JS_weights(c,beta,eps)

alpha = c ./ ( beta + eps ).^2;

sum_alpha = sum(alpha);

w = alpha ./ sum_alpha;

end